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We report the analysis of the statistics of the phase fluctuations in the coda of earthquakes 
recorded during a temporary experiment deployed at Pinyon Flats Observatory, California. The ob- 
served distributions of the first, second and third derivatives of the phase in the seismic coda exhibit 
universal power-law decays whose exponents agree accurately with circular Gaussian statistics. The 
correlation function of the spatial phase derivative is measured and used to estimate the mean free 
path of Rayleigh waves. 
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In the short-period band (> 1 Hz) , ballistic arrivals of 
seismic waves are often masked by scattered waves due to 
small-scale heterogeneities in the lithosphere. The scat- 
tered clastic waves form the pronounced tail of scismo- 
grams known as the seismic coda [T], Q . Even when scat- 
tering is prominent, it is still possible to define the phase 
of the seismic record by introducing the complex analytic 
signal ^{t,r) = j4(t, r)e"^^*''"\ with A the amplitude and 
(f) the phase. In the past, many studies have focused on 
the modeling of the mean field intensity I{t) — {A{t)'^) 
[see H, for review] . The goal of the present paper is to 
study the statistics of the phase field in the coda. In the 
coda, the measured displacements result from the super- 
position of many partial waves which have propagated 
along different paths between the source and the receiver. 
Each path consists of a sequence of scattering events that 
affect the phase of the corresponding partial wave in a 
random way. For narrow-band signals, the phase field can 
therefore be written as (f>{t, r) = ujt + 5(f>{t, r), where uj is 
the central frequency, and S(f) denotes the random fiuctu- 
ations. The trivial cyclic phase tot cancels when a spatial 
phase difference is considered between two neighbouring 
points. Spatially resolved measurements are facilitated 
by dense arrays of seismometers that have been set up 
occasionally. We note that the phase of coda waves has 
not been given much attention so far. The advantage of 
phase is that it is not affected by the earthquake magni- 
tude, and that it contains pure information on scattering, 
not blurred by absorption effects. For the statistical anal- 
ysis of amplitude and phase fluctuations of direct arrivals, 
we refer the reader to e.g. Zheng and Wu 



I. PHASE DISTRIBUTIONS 

We study data sets from a temporary experiment de- 
ployed at Pinyon Flat Observatory (PFO), California, 
in 1990 by an IRIS program. This site exhibits a high 
level of regional seismic activity. The array (see Fig. [T]) 
contained 58 3-components L-22 sensors (2Hz) and was 
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FIG. 1: Geometry of the seismic array 



configured as a grid and two orthogonal arms with sensor 
spacirigs of 7 meters within the grid and 21 meters on the 
arms [5| . We selected 8 earthquakes of magnitude greater 
than 2 with good signal to noise ratio in the coda. Typ- 
ically, epicentral distances are less than 110 km and the 
coda lasts more than 30 seconds after the direct arrivals. 

To perform the statistical analysis, we filtered the sig- 
nal in a narrow frequency band centered around 7 Hz 
(±5%) and selected a 15 s time window starting around 5 
seconds after the direct arrivals. In this time window, the 
signal is believed to be dominated by multiple scattering 
and is highly coherent along the array @. We evaluate 
the Hilbert transform of the vertical displacement which 
yields the imaginary part of the complex analytic signal 
^{t, r) — A{t, r)e"^'^*''"). From the complex field, two defi- 
nitions of the phase can be given: (1) The wrapped phase 
(j) is defined as the argument of the complex field ip in the 
range (— tt : tt]. (2) The unwrapped phase 4>u is obtained 
by correcting for the 27r jumps - occurring when (j) goes 
through the value ±7r - to obtain a continuous function 
with values in R. The (f> distribution is fiat Q. How- 
ever, more information can be extracted by considering 
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higher-order statistics of the phase. For this purpose we 
consider the spatial derivative of the phase, which can be 
estimated in two different ways. 

(1) The first measurement rehes on the difference of 
the wrapped phases A(/) between two seismometers sepa- 
rate by a distance S. Applying the simple finite difference 
formula 0' w A(j)/S an estimate of the spatial derivative 
is obtained. Note that the phase difference At/) takes val- 
ues between —2tt and +27r which does not allow a precise 
estimate for the distribution of the derivative for values 
roughly larger than tt/S. Beyond this value our mea- 
surements will be dominated by finite difference artifacts 
and the distribution is biased by the 27r jumps occurring 
within the distance S. 

(2) The second method uses the difference of the phases 
spatially unwrapped at each time step. This yields 

another estimate of the derivative: (j)' « A(j)u/S which is 
expected to suppress finite difference artifact. In prac- 
tice it is impossible to discriminate a rare but physical 
large phase jump from a small fluctuation that causes a 
27r jump just within the range S. The only possibility 
along 1-D arrays is to impose that the largest admissible 
phase difference between two stations be smaller than tt. 
Hence this (j)' estimate takes values in {—tt/S : tt/5] and is 
biased close to tt/S by the unwrapping processing errors. 
In the limit (5 — > 0, the two definitions ought to be equiva- 
lent because the probability of phase jumps between the 
two stations tends to 0. By averaging over the 8 seis- 
mic records, the lag-time in the coda, the east-west and 
north-south orientations, and the sensor positions within 
the array's grid at fixed S = 7m, we calculate the two re- 
sulting phase derivative distributions which are shown to 
be non-uniform in Figure [2] It is also instructive to con- 
sider the second (third) derivatives of the phase which are 
governed by the 3 (4)-point statistics which are plotted 
in Figure [3l Higher-order derivatives are obtained by ap- 
plying standard finite difference formulas to the wrapped 
phase (this choice is explained in the next section). Since 
the three first derivative distributions are even functions, 
we only represent the positive values. They have all 
similar properties. For small values of the random vari- 
ables, the distributions are nearly flat. For larger values, 
the distributions are governed by a power-law decay ex- 
cept for some peaks which stem from the finite distance 
between the seismic stations. In the following section, 
we will demonstrate that the transition between the two 
behaviours is governed by the wavelength and that the 
power-law decay is a very accurate signature of the Gaus- 
sian nature of the vertical displacements. 

Seismic coda is believed to be composed of multi- 
ply scattered waves. Upon scattering, the many partial 
waves - associated with different paths in the medium 
- would achieve random and independent phase shifts. 
The Gaussian nature would then follow from the cen- 
tral limit theorem. We will thus assert that the coda 
waves obey circular Gaussian statistics (GGS) according 
to which the joint probability of N complex field displace- 
ments ipi, recorded at positions r^, is written as 
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FIG. 2: Distribution of the first derivative of the phase nor- 
malized by the inter-station distance S and measured using 
finite difference formula, (o): wrapped phase; (+): un- 
wrapped phase. The colour lines represent the fits with 
Gaussian theory. Green: wrapped phase; blue: unwrapped 
phase; red : phase derivative. The fitting parameters are 
Q = 2.774 10-" and g{5) = 0.993204. 
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— ii^i^j) is the covariance matrix 
venient to use normalized fields so that C,, 



It is con- 
1. Then, 

the off-diagonal elements are equal to the field correlation 
function = C{\ ri — rj |). From the joint distribution 
of two fields, the probability distribution of the phase 
difference at two points located S apart can be obtained 
by integrating over the amplitudes [9|: 
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where F = gcos{A4>) and g ~ {tpli — 5/2)ip * [r + 5/2)); 
P ^ (27r - |A0|)/47r2 if Acj) is the difference of the 
wrapped phase and P = 1/2-k if A</) is the difference 
of the unwrapped phase Acfi^ ■ In the limit 5 ^ oo of to- 
tally uncorrclatcd fields P(A0) = P. In the limit i5 — > 0, 
we get the following formula for the phase derivative [l^l : 
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where Q = -C"(0) = 2(1 - g)/S^ for 6^0. Figure 
[2] demonstrates that the phase difference between two 
nearby seismometers follows the prediction of Gaussian 
circular statistics with excellent accuracy: we observe a 
good agreement between the theoretical distribution of 
(f)' (Eq. [3]) and the measurements over 3 orders of magni- 
tude in probability and 2 orders of magnitude in deriva- 
tive. A clear discrepancy occurs for large values - typi- 
cally when (j)' tt/S - which can be perfectly explained 
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by the finite distance between the seismic stations. This 
is demonstrated in the same plot which shows excellent 
agreement between formula ([2]) and the measured finite- 
difference statistics of the phase. We observe that the 
formula for the derivative ([3]) agrees with the observa- 
tions on a larger range when the estimate of </>' is based 
on the wrapped phase difference. As a consequence, we 
prefer this method to find the higher derivatives. 

In the frequency band of interest, there is experimen- 
tal evidence that the vertical component of the coda is 
dominated by scattered Rayleigh surface waves [1]. As 
a consequence, we expect the correlation function of the 
field to be given by [lH 

C{r) = (V'(O)V'* (r)) = Jo(fcr) exp(-r/2^), (4) 

which agrees well with observations [3]. Equation [4] con- 
tains two length scales: the wavelength 27r/fc and the 
scattering mean free path (.. We emphasize that formula 
(HI) is valid even when the medium is slightly inelastic. 
Weak absorption does not enter the correlation function 
obtained from the long-time-tail of a diffuse wavefield 
. The form of the correlation function ^ implies that 
Q ~ if fc^ ^ 1. Using the parameter Q obtained by 

fitting the data with equation ([3]), we infer a dominant 
wavelength A of the order of 267 m which agrees with the 
one predicted on the basis of the vertical profile of the 
elastic constants below the PFO array Q offers an 
accurate way of estimating the wavelength, alternative 
to SPAC measurements The use of a narrow band 
signal is crucial because the parameters g and Q strongly 
depend on frequency. 

From the joint Gaussian distribution of 3 and 4 fields, 
we can derive analytically the joint probability functions 
P{4>' 1 4>"Ai P{4>' ! 4'" 1 4'"') I featuring two new constants R 
and S |9|. From these formulas, the marginal distribu- 
tions P{(t)"), P{4)"') can be evaluated numerically. The 
probability distributions of the first, second, and third 
derivatives of the phase exhibit an asymptotic power-law 
decay with exponents —3, —2, —5/3, respectively. These 
universal exponents (i.e. independent of the medium 
properties) can be obtained analytically and provide a 
sensitive fit-independent test of the Gaussian nature of 
the seismic coda. This is illustrated in Figure [3] The 
agreement leaves no doubt that coda waves are in the 
multiple scattering regime. 

II. PHASE DIFFERENCE CORRELATIONS 

We have shown that the distributions of the phase 
derivatives provide accurate information on the short- 
range correlation properties of the field. In the follow- 
ing we use the phase difference correlations to put some 
constraints on the degree of heterogeneity which is re- 
sponsible for long-range correlations along the array. For 
Gaussian statistics and for surface waves obeying Eq. ([5]), 
the phase derivative correlation function is |12l | : 




FIG. 3; Comparison of the Gaussian theory and observations 
for the phase derivatives distribution P(0''"') where n = 1 or 
2 denotes the n*'' derivative of with respect to the spatial 
coordinate. From the fit we find: Q = 2.774 10"" R = 

2.75 10"^ m"^ and S = 4.0 10"^ m"^ 



C^,{r) = i^'m'ir)) = [k/nr) exp(-r/£) (5) 

for r > A. Formula ([5|) has one crucial property. Con- 
trary to the field correlation function C, C^' does not 
oscillate on the wavelength scale but decreases with the 
mean free path as the sole characteristic length scale. 
Any determination of the mean free path based on for- 
mula ([¥]) is impossible because the exponential decay is 
masked by the rapid cyclic oscillations. 

As above, we estimate the phase derivative correlation 
function in a finite difference approximation using the 
formula C^^r) ~ ((A(/)„(r')A0„(r")>/<52 = CA0„(|r' - 
r"|)/5^. Contrary to the probability distribution, we 
found out that the unwrapped phase difference offers a 
better estimate for the phase derivative. The wrapped 
phase difference correlation function is too much domi- 
nated by large 27r jumps with no physical interest. 

The unwrapped phase difference correlation Ca0„ {r) is 
measured along the two orthogonal arms with an aper- 
ture of 252 m and J = 21 m interstation distance. The 
data are averaged over orientation, lag-time in the coda 
and seismic events. The result is presented in linear scale 
in the main part of Figure 3] and shows a decay domi- 
nated by the 1/r factor along the arms of the array, as 
predicted by formula ([5]). This supports our 2D picture 
of wave diffusion. Due to the finite difference C^i (r = 0) 
achieves a finite value at r = 0, which we found to be con- 
sistent with the variance of the unwrapped phase differ- 
ence calculated from (A(/)J) = /^^ (iA^„ (A^„)^P(A</)„) 
and g{5) = 0.98. The parameter g{5) has been deter- 
mined independently by fitting the observed distribution 
of P(A0„) with formula ([2]) for (5 = 21 m. 

Different reasons exist why it is more difficult to mea- 
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FIG. 4: Unwrapped phase difference correlation function, (o): 
experimental results. Dashed line: 1/r fit. Inset: logarithm 
of the correlation multiplied by r for data (o with error bars) 
and numerical simulations at fixed k and £ = 10 km (green), 
£=lkm (red), ^ = 500 m (red) and £ = 200 m (yellow). 



sure the correlation function of the phase derivative than 
e.g. the probability distributions. First, 4th-order statis- 
tics require more averaging to suppress unwanted fluc- 
tuations in the data. Secondly, the intcrstation dis- 
tance 5 = 21 m along the arms reduces the correlation 
between the fields at two nearby stations significantly, 
which favours systematic errors in the derivative. Fi- 
nally, due to frequent breakdowns of the sensors located 
near the ends of the arms, the data could not be aver- 
aged over all sensor positions and all seismic events. As 
a consequence the correlations for distances r > 180 m 
had to be excluded. 

Since formula [5] is valid only in the limit 5 ^ 0, we 
have evaluated the impact of the finite inter-station 
distance inherent to our experimental set-up. To this 
end, we have simulated N correlated Gaussian random 
field displacements on a virtual array with the PFO 



geometry. The results for different values of i at fixed 
k are shown in Figure [4] together with the experimental 
results. The simulated function log(r x Ca</,„) exhibits 
a linear decay with a slope —l/£ as was seen for 5 ^ 
in Eq. (O (see inset in Figure [4]). Therefore, correlations 
of finite-size phase difference still offer direct access to 
the mean free path of the waves in the crust. It is quite 
encouraging to see (inset Figure 2]) that the asymptotic 
exponential regime is already reached for r > A/5. We 
find that the —l/£ slope could in principle be measured 
if the aperture of the network would have been a few 
wavelengths in size (> 500 m), which in general is much 
smaller than the mean free path. Unfortunately, the 
experimental error bars of our data set are too large to 
permit accurate estimates of £ at PFO although it can 
roughly be bounded between I km and 10 km, much 
smaller than the typical path (40 km) taken by coda 
waves coda during 20 s. This has to be interpreted as a 
rough estimate for the mean free path of the Rayleigh 
waves. 



In conclusion, from their observed first three spatial 
derivatives of phase, seismic coda waves are proved to 
obey Gaussian statistics with high accuracy, with a lo- 
cal correlation on the scale of the wavelength. At longer 
scales, we demonstrate that the correlation function of 
the spatial derivative of phase offers a new, promising 
opportunity to measure directly the mean free path £ of 
seismic waves. This measure is not affected by absorp- 
tion neither to scattering anisotropy. As opposed to the 
coherent backscattering effect [lj| , the proposed method 
does not require the sensors to be located in the near 
field of the source, but requires an array of seismic sta- 
tions with sub-wavelength spacing between stations and 
with a total aperture of a few wavelengths. Our work 
highlights the phase as a useful physical object to study 



seismic coda waves. 
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